%% Complex Numbers

i

C=1+i

real(C)

imag(C)

angle(C)

abs(C)

abs(C)*exp(i*angle(C))

%% Waves and images
Xvec=(0:.1:10);

lambda=1;
k=2*pi/lambda;

%Express wave as exponential
Yvec=exp(i*k*Xvec);

%plot the real part
figure
plot(Xvec,real(Yvec))
xlabel('X')
ylabel('Y')

%Show a plane wave in an image
Xmax=25
Ymax=25
[X,Y]=meshgrid((0:.1:Xmax),(0:.1:Ymax))

%look at our X and Y grids
figure
imshow(X,[])
figure
imshow(Y,[])

%Our wave as an image
Wave=exp(i*k*X);
imshow(real(Wave),[])

%Use a wave vector

lambda=2;
Kvec=[2*pi/lambda, 2*pi/lambda]

Xmax=25
Ymax=25
[X,Y]=meshgrid((0:.1:Xmax),(0:.1:Ymax));

Wave=exp(i*(Kvec(1)*X+Kvec(2)*Y));
figure
ImH=imshow(real(Wave),[])
title('e^{\i vec{k}\dot}','interpreter','latex')

%Adding two waves

lambda=2;
Kvec=[2*pi/lambda, 2*pi/lambda]

Xmax=25
Ymax=25
[X,Y]=meshgrid((0:.1:Xmax),(0:.1:Ymax));

Wave=exp(i*(Kvec(1)*X+Kvec(2)*Y));
figure
ImH=imshow(real(Wave),[])

lambda=2;
Kvec2=[2*pi/lambda, -2*pi/lambda]

Xmax=25
Ymax=25
Wave2=exp(i*(Kvec2(1)*X+Kvec2(2)*Y));
figure
ImH=imshow(real(Wave2),[])

Waves=Wave+Wave2;
figure
ImH=imshow(real(Waves),[])

w=5;

h=figure
for tt=1:.1:1000
    figure(h)
    Wave=exp(i*((Kvec(1)*X+Kvec(2)*Y)-w*tt) );
    imshow(real(Wave),[])
    drawnow
    pause(.1)  
end

%% Two Point Sources

lambda=1;
d=20;
w=10;

Xmax=200
Ymax=50

[X1,Y1]=meshgrid((0:.1:Xmax),(-Ymax-d/2:.1:Ymax-d/2));
[X2,Y2]=meshgrid((0:.1:Xmax),(-Ymax+d/2:.1:Ymax+d/2));

k=2*pi/lambda;
R1=sqrt(X1.^2+Y1.^2);
R2=sqrt(X2.^2+Y2.^2);

h=figure
for tt=1:.1:1000
    figure(h)
    Wave1=exp(i*(k*R1-w*tt));% ./R1; %From Huygens Principle
    Wave2=exp(i*(k*R2-w*tt));% ./R2;
    Waves=real(Wave1+Wave2);
    imshow(real(Waves),[])
    drawnow
    pause(.01)  
end


%% Time Average Intensity

Tvec=linspace(0,1/(2*pi*w),20);
WaveAmp=zeros(size(R1));
for tt=1:length(Tvec)
    t=Tvec(tt);
    Wave1=exp(i*(k*R1-w*t));% ./R1; %From Huygens Principle
    Wave2=exp(i*(k*R2-w*t));% ./R2;
    Waves=Wave1+Wave2;
    WaveAmp=WaveAmp+Waves;
end

WaveAmp=WaveAmp/length(Tvec);
figure
imshow((WaveAmp),[])

WaveInt=abs(WaveAmp).^2;
figure
imshow(WaveInt,[])

[X,Y]=meshgrid((0:.1:Xmax),(-Ymax:.1:Ymax));
YScreen=Y(:,end)
IntScreen=WaveInt(:,end)
figure
plot(YScreen,IntScreen)

%From grating equation, d*sin(theta)=m*lambda

L=Xmax;
hold on
for mm=-4:4
    Yfringe=mm*lambda*L/d
    plot([Yfringe,Yfringe],[0 max(IntScreen)],'k--') 
end















